//Name:Chen Xu
//ID:X287284
//Email:xu.c@icloud.com
//Ray tracing global render,use super sampling
//Implement features: 3 different surface(diffuse, mirror reflection, glass)



#include <math.h>
#include <stdlib.h>
#include <stdio.h>

struct Vec {
    double x, y, z;                  // position, also color (r,g,b)
    Vec(double x_=0, double y_=0, double z_=0){ x=x_; y=y_; z=z_; }
    Vec operator+(const Vec &b) const { return Vec(x+b.x,y+b.y,z+b.z); }
    Vec operator-(const Vec &b) const { return Vec(x-b.x,y-b.y,z-b.z); }
    Vec operator*(double b) const { return Vec(x*b,y*b,z*b); }
    Vec mult(const Vec &b) const { return Vec(x*b.x,y*b.y,z*b.z); }
    Vec& norm(){ return *this = *this * (1/sqrt(x*x+y*y+z*z)); }
    double length(){return sqrt(x*x+y*y+z*z);}
    double dot(const Vec &b) const { return x*b.x+y*b.y+z*b.z; } // cross:
    Vec operator%(Vec&b){return Vec(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);}
};

struct Ray {
    Vec o, d;//pos,direction
    Ray(Vec o_, Vec d_) : o(o_), d(d_) {}
};

enum Refl_t {
    DIFF, //solid
    SPEC, //reflect
    REFR  //glass
};  // material types, used in radiance()

struct Sphere {
    double rad;       // radius
    Vec p, e, c;      // position, emission, color
    Refl_t refl;      // reflection type (DIFFuse, SPECular, REFRactive)
    
    //constructer
    Sphere(double rad_, Vec p_, Vec e_, Vec c_, Refl_t refl_):
    rad(rad_), p(p_), e(e_), c(c_), refl(refl_) {}
    
    //check if collsion
    double intersect(const Ray &r) const { // returns distance, 0 if nohit
        // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0
        Vec op = p-r.o;
        double t, eps=1e-4, b=op.dot(r.d), det=b*b-op.dot(op)+rad*rad;
        
        //ray miss sophere
        if (det<0)
            return 0;
        else
            det=sqrt(det);
        return (t=b-det)>eps ? t : ((t=b+det)>eps ? t : 0);
    }
};

struct Plane {
    double w;  //width
    Vec center;
    Vec headup;
    Vec normal;
    Vec color;
    Refl_t refl;
    Plane(double w_,Vec center_,Vec headup_,Vec normal_,Vec color_,Refl_t refl_):
    w(w_),center(center_),headup(headup_),normal(normal_),color(color_),refl(refl_){}
    
};

//  Draw scene,every object is based on sphere
//  Scene: radius, position, emission, color, material
Sphere spheres[] = {//Scene: radius, position, emission, color, material
    //  Sphere(1e5, Vec( 1e5+1,40.8,81.6), Vec(),Vec(.75,.25,.25),SPEC),//Left
    //  Sphere(1e5, Vec(-1e5+99,40.8,81.6),Vec(),Vec(.25,.25,.75),SPEC),//Rght
    // Sphere(1e5, Vec(50,40.8, 1e5),     Vec(),Vec(.75,.75,.75),DIFF),//Back
    // Sphere(1e5, Vec(50,40.8,-1e5+170), Vec(),Vec(),           DIFF),//Frnt
    Sphere(1e5, Vec(50, 1e5, 81.6),    Vec(),Vec(.75,.75,.75),DIFF),//Botm
    // Sphere(1e5, Vec(50,-1e5+81.6,81.6),Vec(),Vec(.75,.75,.75),DIFF),//Top
    
    Sphere(16.5,Vec(27,16.5,47),       Vec(),Vec(1,1,1)*.999, SPEC),//Mirr
    Sphere(16.5,Vec(73,16.5,78),       Vec(),Vec(1,1,1)*.999, DIFF),//Glas
    
    Sphere(1.5, Vec(50,81.6-16.5,81.6),Vec(4,4,4)*100,  Vec(), DIFF),//Lite
};
//bound x to 0-1
inline double clamp(double x){ return x<0 ? 0 :(x>1 ? 1 : x); }

//convert color to 0-1 to range 0-255
inline int toInt(double x){ return int(pow(clamp(x),1/2.2)*255+.5); }


//check whether this sphere intersect with ray r
inline bool intersect(const Ray &r, double &t, int &id){
    
    //number of sphreres
    double n=sizeof(spheres)/sizeof(Sphere);
    double d;
    double inf=t=1e20;
    
    for(int i=int(n);i--;){
        if((d=spheres[i].intersect(r)) && d<t){
            t=d;
            id=i;
        }
    }
    return t<inf;
}
int numSpheres = sizeof(spheres)/sizeof(Sphere);

Vec radiance(const Ray &r, int depth, unsigned short *Xi){
    double t;                               // distance to intersection
    int id=0;                               // id of intersected object
    if (!intersect(r, t, id))
        return Vec(); // if miss, return black
    
    const Sphere &obj = spheres[id];        // the hit object
    
    Vec x=r.o+r.d*t;                        //point of collison,head of light,move distance t
    Vec n=(x-obj.p).norm();                 //normal vector
    Vec nl=n.dot(r.d)<0?n:n*-1;             //钝角或锐角->n正反
    Vec f=obj.c;                            //color
    
    double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z; // max refl
    
    
    //hard code recursion depth
    if (++depth>8) return obj.e;
    
    // Ideal DIFFUSE reflection
    if (obj.refl == DIFF){
        double r1=M_PI_2*erand48(Xi);  // random angle
        double r2=erand48(Xi);         // random distance
        double r2s=sqrt(r2);
        Vec w = nl;                     //normalized
        Vec u = ((fabs(w.x)>.1?Vec(0,1):Vec(1))%w).norm();
        Vec v = w%u;
        Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm();
        return obj.e + f.mult(radiance(Ray(x,d),depth,Xi));
    }
    // Ideal mirror reflection
    else if (obj.refl == SPEC)
        return obj.e + f.mult(radiance(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi));
    
    // Ideal dielectric REFRACTION
    else{
        //compute reflect ray
        Ray reflRay(x, r.d-n*2*n.dot(r.d));
        bool into = n.dot(nl)>0;                // Ray from outside going in?
        double nc=1, nt=1.5;
        double nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
        if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0)    // Total internal reflection
            return obj.e + f.mult(radiance(reflRay,depth,Xi));
        Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();
        double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n));
        double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P);
        return obj.e + f.mult(depth>2 ? (erand48(Xi)<P ? radiance(reflRay,depth,Xi)*RP:radiance(Ray(x,tdir),depth,Xi)*TP) : radiance(reflRay,depth,Xi)*Re+radiance(Ray(x,tdir),depth,Xi)*Tr);
    }
}

Vec radiance1(const Ray &r_, int depth_, unsigned short *Xi){
    double t;                               // distance to intersection
    int id=0;                               // id of intersected object
    Ray r=r_;
    int depth=depth_;
    // L0 = Le0 + f0*(L1)
    //    = Le0 + f0*(Le1 + f1*L2)
    //    = Le0 + f0*(Le1 + f1*(Le2 + f2*(L3))
    //    = Le0 + f0*(Le1 + f1*(Le2 + f2*(Le3 + f3*(L4)))
    //    = ...
    //    = Le0 + f0*Le1 + f0*f1*Le2 + f0*f1*f2*Le3 + f0*f1*f2*f3*Le4 + ...
    //
    // So:
    // F = 1
    // while (1){
    //   L += F*Lei
    //   F *= fi
    // }
    Vec cl(0,0,0);   // accumulated color
    Vec cf(1,1,1);  // accumulated reflectance
    while (1){
        if (!intersect(r, t, id)) return cl; // if miss, return black
        const Sphere &obj = spheres[id];        // the hit object
        Vec x=r.o+r.d*t, n=(x-obj.p).norm(), nl=n.dot(r.d)<0?n:n*-1, f=obj.c;
        double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z; // max refl
        cl = cl + cf.mult(obj.e);
        if (++depth>5) if (erand48(Xi)<p) f=f*(1/p); else return cl; //R.R.
        cf = cf.mult(f);
        if (obj.refl == DIFF){                  // Ideal DIFFUSE reflection
            double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), r2s=sqrt(r2);
            Vec w=nl, u=((fabs(w.x)>.1?Vec(0,1):Vec(1))%w).norm(), v=w%u;
            Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm();
            //return obj.e + f.mult(radiance(Ray(x,d),depth,Xi));
            r = Ray(x,d);
            continue;
        } else if (obj.refl == SPEC){           // Ideal SPECULAR reflection
            //return obj.e + f.mult(radiance(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi));
            r = Ray(x,r.d-n*2*n.dot(r.d));
            continue;
        }
        Ray reflRay(x, r.d-n*2*n.dot(r.d));     // Ideal dielectric REFRACTION
        bool into = n.dot(nl)>0;                // Ray from outside going in?
        double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
        if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0){    // Total internal reflection
            //return obj.e + f.mult(radiance(reflRay,depth,Xi));
            r = reflRay;
            continue;
        }
        Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();
        double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n));
        double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P);
        // return obj.e + f.mult(erand48(Xi)<P ?
        //                       radiance(reflRay,    depth,Xi)*RP:
        //                       radiance(Ray(x,tdir),depth,Xi)*TP);
        if (erand48(Xi)<P){
            cf = cf*RP;
            r = reflRay;
        } else {
            cf = cf*TP;
            r = Ray(x,tdir);
        }
        continue;
    }
}
Vec radiance2(const Ray &r, int depth, unsigned short *Xi,int E=1){
    double t;                               // distance to intersection
    int id=0;                               // id of intersected object
    if (!intersect(r, t, id)) return Vec(); // if miss, return black
    const Sphere &obj = spheres[id];        // the hit object
    Vec x=r.o+r.d*t, n=(x-obj.p).norm(), nl=n.dot(r.d)<0?n:n*-1, f=obj.c;
    double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z; // max refl
    if (++depth>5||!p) if (erand48(Xi)<p) f=f*(1/p); else return obj.e*E;
    if (obj.refl == DIFF){                  // Ideal DIFFUSE reflection
        double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), r2s=sqrt(r2);
        Vec w=nl, u=((fabs(w.x)>.1?Vec(0,1):Vec(1))%w).norm(), v=w%u;
        Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm();
        
        // Loop over any lights
        Vec e;
        for (int i=0; i<numSpheres; i++){
            const Sphere &s = spheres[i];
            if (s.e.x<=0 && s.e.y<=0 && s.e.z<=0) continue; // skip non-lights
            
            Vec sw=s.p-x, su=((fabs(sw.x)>.1?Vec(0,1):Vec(1))%sw).norm(), sv=sw%su;
            double cos_a_max = sqrt(1-s.rad*s.rad/(x-s.p).dot(x-s.p));
            double eps1 = erand48(Xi), eps2 = erand48(Xi);
            double cos_a = 1-eps1+eps1*cos_a_max;
            double sin_a = sqrt(1-cos_a*cos_a);
            double phi = 2*M_PI*eps2;
            Vec l = su*cos(phi)*sin_a + sv*sin(phi)*sin_a + sw*cos_a;
            l.norm();
            if (intersect(Ray(x,l), t, id) && id==i){  // shadow ray
                double omega = 2*M_PI*(1-cos_a_max);
                e = e + f.mult(s.e*l.dot(nl)*omega)*M_1_PI;  // 1/pi for brdf
            }
        }
        
        return obj.e*E+e+f.mult(radiance2(Ray(x,d),depth,Xi,0));
    } else if (obj.refl == SPEC)              // Ideal SPECULAR reflection
        return obj.e + f.mult(radiance2(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi));
    Ray reflRay(x, r.d-n*2*n.dot(r.d));     // Ideal dielectric REFRACTION
    bool into = n.dot(nl)>0;                // Ray from outside going in?
    double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
    if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0)    // Total internal reflection
        return obj.e + f.mult(radiance2(reflRay,depth,Xi));
    Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();
    double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n));
    double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P);
    return obj.e + f.mult(depth>2 ? (erand48(Xi)<P ?   // Russian roulette
                                     radiance2(reflRay,depth,Xi)*RP:radiance2(Ray(x,tdir),depth,Xi)*TP) :
                          radiance2(reflRay,depth,Xi)*Re+radiance2(Ray(x,tdir),depth,Xi)*Tr);
}

Vec radiance3(const Ray &r, int depth, unsigned short *Xi){
    double t;                               // distance to intersection
    int id=0;                               // id of intersected object
    if (!intersect(r, t, id)) return Vec(); // if miss, return black
    const Sphere &obj = spheres[id];        // the hit object
    Vec x=r.o+r.d*t, n=(x-obj.p).norm(), nl=n.dot(r.d)<0?n:n*-1, f=obj.c;
    double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z; // max refl
    if (++depth>5) if (erand48(Xi)<p) f=f*(1/p); else return obj.e; //R.R.
    if (obj.refl == DIFF){                  // Ideal DIFFUSE reflection
        double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), r2s=sqrt(r2);
        Vec w=nl, u=((fabs(w.x)>.1?Vec(0,1):Vec(1))%w).norm(), v=w%u;
        Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm();
        return obj.e + f.mult(radiance3(Ray(x,d),depth,Xi));
    } else if (obj.refl == SPEC)            // Ideal SPECULAR reflection
        return obj.e + f.mult(radiance3(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi));
    Ray reflRay(x, r.d-n*2*n.dot(r.d));     // Ideal dielectric REFRACTION
    bool into = n.dot(nl)>0;                // Ray from outside going in?
    double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
    if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0)    // Total internal reflection
        return obj.e + f.mult(radiance3(reflRay,depth,Xi));
    Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();
    double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n));
    double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P);
    return obj.e + f.mult(depth>2 ? (erand48(Xi)<P ?   // Russian roulette
                                     radiance3(reflRay,depth,Xi)*RP:radiance3(Ray(x,tdir),depth,Xi)*TP) :
                          radiance3(reflRay,depth,Xi)*Re+radiance3(Ray(x,tdir),depth,Xi)*Tr);
}


void mymethod(){
    int w=640, h=640;
    int samplewidth = 5;
    int samps = samplewidth*samplewidth;                        // # samples
    Vec cam(50,50,120);
    Vec screen(48,48,119);
    double factor = 4/w;
    
    //Ray cam(Vec(50,52,295.6), Vec(0,-0.042612,-1).norm());      // camera pos, dir
    //Vec cx=Vec(w*.5135/h);                                      // x direction increment
    //Vec cy=(cx%cam.d).norm()*.5135;                             // y direction increment
    
    //Vec r;                                                      // colors for samples
    Vec *c=new Vec[w*h];                                        // store the image
    
    for (int y=0; y<h; y++){                                    // Loop over image rows
        //   fprintf(stderr,"\rRendering (%d spp) %5.2f%%",samps*4,100.*y/(h-1));
        for (unsigned short x=0, Xi[3]={0,0,static_cast<unsigned short>(y*y*y)}; x<w; x++)  { // Loop cols
            for (int jj= 0; jj<samplewidth; jj++) {
                Vec dirction = Vec(x-w/2,y-h/2,-400);
                Ray ray = Ray(cam+dirction.norm(),dirction.norm());
                //  fprintf(stdout, "d:%f,%f,%f\n",ray.d.x,ray.d.y,ray.d.z);
                Vec r = Vec();
                //accumulate color
                r = r + radiance(ray,0,Xi);
                int i = y*w+x;
                // fprintf(stdout, "d:%f,%f,%f\n",d.x,d.y,d.z);
                
                //fprintf(stdout, "r:%f,%f,%f\n",r.x,r.y,r.z);
                c[i] = c[i] + Vec(clamp(r.x),clamp(r.y),clamp(r.z))*.25;
                
            }
            
        }
    }
    FILE *f = fopen("/Users/chenxu/desktop/myppm.ppm", "w");         // Write image to PPM file.
    fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);
    for (int i=0; i<w*h; i++)
        fprintf(f,"%d %d %d ", toInt(c[i].x), toInt(c[i].y), toInt(c[i].z));
}

void example(){
    int samp = 3;
    int w=640, h=640, samps =samp*samp; // # samples
    Ray cam(Vec(50,52,295.6), Vec(0,-0.042612,-1).norm()); // cam pos, dir
    Vec cx=Vec(w*.5135/h), cy=(cx%cam.d).norm()*.5135, r, *c=new Vec[w*h];
#pragma omp parallel for schedule(dynamic, 1) private(r)       // OpenMP
    for (int y=0; y<h; y++){                       // Loop over image rows
        fprintf(stderr,"\rRendering (%d spp) %5.2f%%",samps*4,100.*y/(h-1));
        for (unsigned short x=0, Xi[3]={0,0,static_cast<unsigned short>(y*y*y)}; x<w; x++)  // Loop cols
            for (int sy=0, i=(h-y-1)*w+x; sy<samp; sy++)                                    //  subpixel rows
                for (int sx=0; sx<samp; sx++, r=Vec()){                                     //  subpixel cols
                    for (int s=0; s<samps; s++){
                        double
                        r1=2*erand48(Xi), dx=r1<1 ? sqrt(r1)-1: 1-sqrt(2-r1);
                        double r2=2*erand48(Xi), dy=r2<1 ? sqrt(r2)-1: 1-sqrt(2-r2);
                        Vec d = cx*( ( (sx+.5 + dx)/2 + x)/w - .5) +
                        cy*( ( (sy+.5 + dy)/2 + y)/h - .5) + cam.d;
                        r = r + radiance2(Ray(cam.o+d*140,d.norm()),0,Xi)*(1./samps);
                    }
                    c[i] = c[i] + Vec(clamp(r.x),clamp(r.y),clamp(r.z))*.25;
                }
    }
    FILE *f = fopen("/Users/chenxu/desktop/example.ppm", "w");         // Write image to PPM file.
    fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);
    for (int i=0; i<w*h; i++)
        fprintf(f,"%d %d %d ", toInt(c[i].x), toInt(c[i].y), toInt(c[i].z));
}

int main(int argc, char *argv[]){
    example();
}